function [IRF_Gibbs, AA] = IRF_BVARGibbs(gamma_g,Su_g,p,N,Gibbs,horizon);
%%
for jg = 1:Gibbs
    
    A = zeros(N*p);
    A(1:N,:) = gamma_g(1:end-1,:,jg)';
    A(N+1:end,1:N*(p-1)) = eye(N*(p-1));
    J = zeros(N,N*p); J(:,1:N) = eye(N);
    
    B0 = chol(Su_g(:,:,jg))'; %% Lower triangular
    AA(:,:,jg)=A;
    for j = 0:horizon
        IRF_Gibbs(:,:,j+1,jg) = J*A^j*J'*B0;
    end     
end
